A perturbative approach to the Bak-Sneppen Model 
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We study the Bak-Sneppen model in the probabilistic framework of the Run Time Statistics (RTS). 
This model has attracted a large interest for its simplicity being a prototype for the whole class 
of models showing Self-Organized Criticality. The dynamics is characterized by a self-organization 
of almost all the species fitnesses above a non-trivial threshold value, and by a lack of spatial and 
temporal characteristic scales. This results in avalanches of activity power law distributed. In this 
letter we use the RTS approach to compute the value of x c , the value of the avalanche exponent r 
and the asymptotic distribution of minimal fitnesses. 
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The concept of Self-Organized Criticality (SOC) has 
been introduced in order to explain the ubiquitous pres- 
ence of scale free phenomena in Nature. Ranging from 
fractures jjj to river basins Q, the lack of spatial and/or 
temporal characteristic scales is the main feature of these 
processes. To detect if an underlying general mechanism 
exists, two paradigmatic models have been introduced. 
One mimics the evolution of sandpiles || , and the other 
one, has been originally introduced to mimic punctuated 
equilibrium in the evolution of the ecology of species. 
We focus here on the latter one, known as Bak-Sneppen 
(BS) model Q. We develop a perturbative approach to 
BS based on the Run Time Statistics (5|(|. This allows 
us to evaluate the main quantities for the model in the 
stationary state: the fitnesses histogram, the distribution 
of minimal fitnesses, and the avalanche exponent. 

The BS model is defined as follows: an 1-d ring of 
N different sites models the N species present in the 
ecosystem. Links between first nearest neighbor (fnn) 
sites model interaction between species (i.e. predation). 
A random number Xi G [0,1], extracted from the uniform 
probability density function (pdf) fo(x) — 1, is assigned 
to each site/species i in the system. The x^s quantify the 
fitness of species i to survive the competition with the 
other species in the environment. At each time step, the 
species with the minimal fitness is selected and removed. 
The species connected to it through predation (i.e. the 
fnn) are also removed, regardless their fitness. Three 
new species with randomly extracted fitnesses take their 
place. We refer to this updating rule of the model as the 
refreshing rule. After a transient period this system self- 
organizes in a stationary state characterized by two main 
features: (1) the normalized fitnesses histogram (j)(x) is 
given, apart from corrections of order 1/N, by 4>{x) — 
9(x — x c )/(l — x c ), where 9{x) is the usual step function, 
and x c — 0.66702(3) 0; (2) the stationary dynamics 
evolves as a sequence of critical avalanches [^J^,^], the 
duration s of which is power law distributed: -P(s) ~ s~ T 
with t ~ 1.07. 

The best analytical results for BS come from mean field 
approaches giving a threshold x c = 1/3 JTo| , very far from 
the real numerical value. This striking disagreement lead 



us to investigate if the actual value can be computed 
by using a theoretical tool, called Run Time Statistics 
(RTS) |,§. RTS is designed to study growth processes 
in a medium with quenched disorder. In particular, the 
main results jll],[l2j obtained through RTS deal with the 
class of models derived by the Invasion Percolation (IP) 
in d — 2 fLU , to which the BS model can be related. 

RTS is based on the intuitive idea that the larger the 
number of time-steps one species survives, then probably 
the larger its fitness is. In other words, if at the first 
time-step with no information, one assigns a uniform pdf 
for the species fitnesses, at successive time-steps informa- 
tion about the history increases the effective pdf at high 
fitness values for the surviving species. More formally, 
this information is stored in effective (conditional) pdf 's 
{fi,t(x)} of the quenched disorder variables {x^}. These 
effective pdf 's can then be used to develop a step by step 
algorithm to describe probabilistically the whole tree of 
the possible dynamics of the system, starting from the 
last available situation ||. Let us suppose to know at 
time t the set of effective pdf's {fi,t(%)}- By using joint 
probability, it is possible to write for each species the ef- 
fective probability to be selected at that time-step, and 
one can update the effective pdf's, obtaining the new set 
of functions {fi^+i(x)}. For example, consider a limit 
case with a system composed only by two sites a, b. The 
minimum is selected and refreshed without affecting the 
other. At the beginning (t — 0) both pdf's are uniform 
and equal to 1, and the probability /ia.oiMM to be se- 
lected for sites a,b are both equal to 1/2. We assume 
without loss of generality that x a < Xb- By using the 
postulate of the theory of conditioned probability 



fa, i( x a) = 1 

fb, l( X b) = ^ ^ [ dx a fafi{x a ) = 2x b 



Ma,0 Jo 

and consequently 

Mo,X = / dXbfb,l(Xb) / dx a f a ,i(x a ) = 2/3 
/o Jo 

jti M = / dx a f a .i{x a ) / dx b fb,i{xb) = 1/3. 
Jo Jo 
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FIG. 1. the first three temporal steps of an avalanche tree. 
The initiator i — is selected at t = and the possible 
growths are shown. 



That is the just refreshed site has a larger probability to 
grow. This argument can be easily generalized for any 
time t. For instance, if the site selected is always the 
same, / M = (t + and fj,b,t = — T 

In order to implement this program for the BS, one 
needs only an initial site configuration with known 
{fi,t(x)}. There are two possibilities: (i) to follow the 
dynamics from the first time step t = 0, having obviously 
fi,o(x) = fo(x) = 1 for each i; (ii) to start at an arbitrary 
time t for which the set {fi,t(x)} is known. In the follow- 
ing, we consider mainly the case (ii): the system is consid- 
ered already at the critical stationary state, just after the 
first step of a critical avalanche. We call the site selected 
at this first step the initiator of the avalanche. We recall 
that the initiator fitness lies in an interval of order 1 /N 
around x c , and all the other sites have a fitness x > x c . 
Avalanches are defined as geometrically and causally con- 
nected sequences of growths. Since it has been proved 
that each avalanche is independent on the others, the 
probabilistic study of a generic critical avalanche pro- 
vides a complete statistical description of the stationary 
state. Without loss of generality, we consider the initia- 
tor placed at i = and its selection at t — 0. At t = 1, 
the just grown initiator and its two fnn sites have the fit- 
ness uniformly distributed because of the refreshing rule 
(i.e. foA(x) = f-i.i{x) = fi,i(x) = f (x) = 1). Let 
us call A t (active sites) the set of sites refreshed by the 
avalanche up to its t th step. Consequently, A\ is com- 
posed by the three sites — 1,0, +1 . The set A t defines 
a connected segment on the system because of the geo- 
metrical connectivity of an avalanche. One can represent 
the evolution of the system as a branching process (see 
Fig. Q). At each time-step, the growth of each active 
site corresponds to a branching event. The growth of 
"non-active" sites are not included in the branching tree, 
because this implies the end of the avalanche and the 
beginning of a new one. Any finite connected path on 
this branching tree represents a single realization of the 
avalanche. The tree contains all the possible realizations 
of the avalanche with a given initiator. 

The RTS approach to an avalanche in the critical sta- 
tionary state can be formulated as follows. The avalanche 
proceeds, after the selection of the initiator, if at least 



one of the three sites belonging to A\ has fitness below 
x c . In general, at the t th time-step of the avalanche, if 
one knows the effective pdf fi,t(x) for each site i G A t , 
one can write the effective growth probability of site i 
conditioned to the fixed past history: 



dxif iit (xi) Yl 



(3) 



Each of these possible growth events corresponds to a 
branch of the tree leading from the starting configuration 
to a new one characterized by a new set At^i of active 
sites with their new effective pdf 's. The set of sites At+\ 
is related with the "mother" set At. If the just grown 
site i is not an extreme of the segment defined by A t , 
then A t +i = A t . If, instead, i is the left (right) extreme 
of A t , then the left (right) fnn of i must be added to 
A t in order to obtain A t +\. After the selection of the 
site i, the new effective pdf's of the sites belonging to 
A t+ i are obtained from those of A t by applying both the 
refreshing rule of the dynamics and the law of conditional 
probability P(B\A) = P ^^" > to store the information 
about the last growth event: 

[ Jk,t{Xc) X > X c 

h,t+i(x) =/»_i,tfi(x) = f l+ i, t +i{x) = fo(x) = 1 (4) 

where k is any species in the set A t different from i and its 
neighbors, and j is an active site different from k, i. The 
pdf of an active site is always constant above x c because 
in the stationary state no site with x > x c grows, then 
no information is available for x > x c . Through this 
algorithm one can store step by step information on the 
fitnesses from the actual fixed history of the avalanche. 
By iterating Eqs.(^,^), one can write for a time-length At 
the probability of any connected path C(At) belonging to 
the tree as 



A/ 



w(c(At)) = \: 



(5) 



where it is a site growing at time t th of the path C(At). 

Let us now write an equation for the average fitness 
histogram 4>t{x). From the refreshing rule, one has 



Ncj)t +1 (x) = N(j> t (x) - m i>t+ i{x) - /i_i, i+ i(x) 
-fi+i,t+i(x) +3, 



(G) 



where m itt +i(x), fi^i >t +i(x) and f i+lit (x) are the effec- 
tive pdf's that sites + l would have respectively 
after the selection of site i, if their fitnesses were not 
refreshed. Eq.(||) states that three species are removed 
with their pdf's, and three new species with uniform pdf 
enter. Considering the average (...) of Eq.(||) over all the 
possible paths, and imposing the stationarity condition 
4>t+i{x) = 4>t{x) = 4>{x), one has: 
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(mi(x)) + + fi+i(x)\ 



3, 



(7) 



We use Eq.([7J) to estimate x c and the behavior of 4>(x) 
under the simple hypothesis that all the sites not touched 
by the avalanche have a fitness x > x c (which comes from 
the stationarity condition). The functions fi—i,t+i{x) 
and f i+ i.t+i(x) are given by the second line of Eq.(^) 
with k = i — 1 and k = i + 1 respectively. In order to 
perform the average (...), one has to consider all the ways 
in which the selection of site i represents the t step of 
an avalanche (i.e. one has to consider all the finite paths 
of the tree of Fig. pi). Then one can write: 



fi-i(x) + f i+ i(x) 



EtEc(t)W(C(t)) fi- llt+ i(x) +f i+1 , t+1 (x) 



(8) 



C(t) 



where W(C(t)) is the RTS statistical weight of the generic 
path C(t) of length t, in the tree of Fig. [I], given by 

indicates that the functions in- 



Jc(t; 



Eq.(||). Moreover 
side brackets are evaluated by applying step by step the 
RTS algorithm to the path C(t). Note, anyway, that 
if i is the left (right) extreme of At, then the knowl- 
edge of the path does not give any information about 
the sites i — 1 (i + 1) apart that, before refreshment, 
aft— l > x c (a5j+i > x c ). In this case we can approxi- 
mate [fi-l,t+l{x)] C {t) = <t>( x ) ([fi+l,t+l(x)]c(t) = H x ))- 

Distinguishing these last terms in <j){x) in Eq. (|8|) from 
the others, one can write: 

[fi-i{x) + f i+1 {x)) = A{x c )<t>{x) + B{x,x c ), (9) 

where A(x c ) and B(x,x c ) are positive functions, and 
B(x,x c ) — B(x c ,x c ) for x > x c , because of Eq. (^). 

In principle, we should evaluate (mj(x)) applying the 
rule of conditional probability in analogy with Eq.(|j), 
and then performing the average as in Eq. (JsJ) . However, 
this would reduce Eq.(^) to a trivial identity. For this 
reason, in order to evaluate (m%(x)}, we adopt another 
strategy in the optics of the point (i). We impose that 
Xi is the minimal fitness in the whole system at a certain 
time-step T 3> N from the actual beginning of the dy- 
namics. Consequently, we apply the RTS algorithm from 
the first time-step f| of the dynamics up to T (in this 
case we have to consider all the system sites). In the end, 
we consider the average (...) over all the possible histories 
up to time T. Since for T ^> N the system must be in the 
stationary state H, then, after the average, this descrip- 
tion must be consistent with the previous one where the 
system was considered directly in the stationary state. 
In order to simplify the mathematical treatment, we use 
the formula obtained for the generalized (3-BS model pE) , 
and we take the limit f3 — > oo at the end of calculations: 
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FIG. 2. Empty points represent the values of x c obtained 
by the application of the RTS algorithm from n — 2 to n = 7. 
The continuous line represents the fit curve x c (n) = 0.66 — ax 6 
with a = 0.291 ± 0.003 and b = 0.20 ± 0.03. The insert shows 
the behavior of fl(ra) up to n = 7. Assuming Q(n) ~ n~ T+1 
as a good approximation also for small n, one finds r ~ 1.05. 



where x c is defined by the relation e ^ Xc = 
J dx (j){x)e~ l3x with |9 > 1, which gives also a deeper 
meaning to x c 15 1. Note that in this contest the presence 
of a finite threshold x c arises naturally, without external 
assumption. In order to obtain Eq. (|l0|), we used the 
central limit theorem for N 1 M] and we made the 
following approximation flifl : 
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Introducing Eqs. @ and jlfj| ) in Eq. (Q), we obtain: 

3 - B(x, x c ) 



4>(x) 



A(x c ) + N6(x c - x) 



(11) 



(rrii(x)) = N<f){x)6(x c — x) 



(10) 



Note that, because of the behavior of B(x,x c ), <j)(x) is 
constant for x > x c , while for x < x c is of order l/N. 
There is only one value of x c for which </>(x) is normalized 
(i.e. Jq 1 4>{x)dx = 1). This value gives the value x c of the 
BS threshold. 

In principle, one should consider the contribution to 
Eq. (ph coming from any length t in the avalanche tree 
(Fig.pl). By stopping the sums at a value t = n, we obtain 
an n-order approximation. We have evaluated the prob- 
ability of the paths, through RTS, up to n = 7 by using a 
computer program for numerical integration. The results 
for x c , for n — 2 up to n = 7, are reported in Fig.^j: The 
best evaluation is x c (n — 7) ~ 0.465, much better than 
the mean- field result x c = 1/3, but still quite far from 
the above reported numerical value. However, we veri- 
fied that this behavior is compatible with an asymptotic 
value x c (n — > oo) ~ 0.66. To this aim, we made a fit with 
the simplest possible function x c (n) = 0.66 — ax b (Fig.^) 
which is well compatible with the given asymptotic value. 
The fit vaues are a = 0.291 ± 0.003 and b = 0.20 ± 0.03. 
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FIG. 3. The continuous line gives the theoretical station- 
ary distribution m(x) of the minimal fitnesses where n — 7 
and assuming x c (n — > oo) = 0.66. The points represent the 
numerical behavior evaluated in extensive simulations. 



2 + N0(x c - x) 



(12) 



that is 4>(x) — 0(l/N) for x < x c and 4>(x) = 3/2 for 
x > x c . Normalizing (j)(x), one gets x c = 1/3, which is 
the above cited mean- field result. 

In conclusion, this paper presents a kind of perturba- 
tive approach to the BS model, based on the probabilistic 
framework called Run Time Statistics. Through this ap- 
proach, we compute the self-organized threshold x c , the 
avalanche exponent r, and the stationary distribution of 
minimal fitnesses m(x). In principle one could also ob- 
tain in this way the exponent /i characterizing the sites 
covered by an avalanche. However, since we can reach 
at the best the 7 t h time step of growth, this results in a 
very poor statistics for the number of points covered by 
the paths. Authors acknowledge the support of the EU 
Network ERBFMRXCT980183. 



This small value of b is due to the fact that the avalanche 
size distribution P(s) Q is characterized by a small expo- 
nent (r ~ 1.07), henceforth all the sizes s are important 
for statistics. One can use x c (n — > oo) to evaluate both 
the avalanche exponent r and the average minimum dis- 
tribution m(x) = {m,i{x)) . The exponent r can be found 
studying, using x c (n — > oo) = 0.66 in the Eqs. (^) and 
(|), the behavior of Q(t) = J2c(t) w i c i t )) for t ranging 
from 1 to the maximal n possible. Indeed fl(t) is propor- 
tional to the probability that the avalanche lasts at least t 
time-steps. In the scaling regime, il(t) ~ t~ T+1 . Making 
this hypothesis, the result for n = 7 is r(n = 7) ~ 1.05 
(see insert in Fig.^J) , which is in agreement with the 
known numerical value. Clearly, the hypothesis of scal- 
ing regime for n = 7 is a strong hypothesis, but it is 
based on the fact that for Invasion Percolation, which 
has similar dynamical rules, it has been shown [jll| that 
the microscopical dynamical rules are already the scale 
invariant ones. Finally, we can obtain an approxima- 
tion, of the distribution of the minimal fitness m(x) in 
the stationary state considering the Eqs. ( |l(]| ) and (11), 
where the latter one is evaluated using RTS up to n = 7, 
but using x c (n — > oo) in the RTS calculations. Imposing 
normalization, we obtain the function m(x) reported in 
Fig.^. In the same figure this result is compared with 
the known asymptotic numerical distribution of mini- 
mal fitnesses j|. The agreement is again quite good, 
considering the strong approximation in truncating the 
sums of Eq.(|8|) at n — 7. Finally, it is worth to note 
that our method include also the mean-field results for 
BS. We recall that g(J in the mean-field version of the 
model, when a site i is selected, the other two sites re- 
freshed are not the fnn sites, but two randomly chosen 
system sites. This means that in Eq.(^) one has to use 
(fi-i(x)) = (fi+i(x)) = 4>{x). Henceforth, the asymp- 
totic histogram becomes: 
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